/*==============================================================================
案例3：销售预测分析 - 完整流程
================================================================================

业务场景：
某公司销售部门需要预测未来的销售额，以便：
1. 制定销售目标和计划
2. 优化库存管理
3. 合理分配销售资源
4. 评估营销活动效果

数据说明：
使用 Stata 内置的数据集，模拟销售数据
包含产品特征、市场因素、历史销售等信息

分析方法：
1. OLS 回归（基准模型）
2. Lasso 回归（特征选择）
3. Elastic Net（正则化）
4. 模型对比和选择

作者：Stata ML Course
日期：2025-11-03
==============================================================================*/

clear all
set more off
capture log close
log using "output/cases/case03_sales_forecasting.log", replace

*------------------------------------------------------------------------------
* 1. 数据准备
*------------------------------------------------------------------------------

display "=" * 80
display "案例3：销售预测分析"
display "=" * 80
display ""

* 加载数据（使用 auto.dta 模拟销售数据）
sysuse auto, clear

display "原始数据概览："
describe
display "样本量: " _N

*------------------------------------------------------------------------------
* 2. 特征工程（创建销售相关变量）
*------------------------------------------------------------------------------

display ""
display "特征工程：创建销售预测特征"
display "-" * 80

* 2.1 目标变量：销售额
set seed 20251103
gen sales_amount = price * (10 + int(runiform() * 50))
label variable sales_amount "月销售额（千元）"

* 2.2 产品特征
gen product_quality = 100 - mpg  // 质量得分
label variable product_quality "产品质量得分"

gen product_size = weight / 1000
label variable product_size "产品规格"

gen product_performance = displacement / 100
label variable product_performance "产品性能指数"

* 2.3 价格特征
gen unit_price = price / 1000
label variable unit_price "单价（千元）"

gen price_level = 1 if price < 5000
replace price_level = 2 if price >= 5000 & price < 10000
replace price_level = 3 if price >= 10000 & price < .
label define price_lbl 1 "低价" 2 "中价" 3 "高价"
label values price_level price_lbl
label variable price_level "价格档次"

* 2.4 市场因素
gen market_share = runiform() * 30
label variable market_share "市场份额（%）"

gen competitor_count = int(runiform() * 10) + 1
label variable competitor_count "竞争对手数量"

gen brand_strength = (foreign == 0) * 60 + runiform() * 40
label variable brand_strength "品牌力度"

* 2.5 营销投入
gen marketing_spend = price * runiform() * 0.1
label variable marketing_spend "营销费用（千元）"

gen promotion_intensity = runiform() * 100
label variable promotion_intensity "促销力度"

* 2.6 季节性因素
gen season = int(runiform() * 4) + 1
label define season_lbl 1 "春季" 2 "夏季" 3 "秋季" 4 "冬季"
label values season season_lbl
label variable season "销售季节"

* 2.7 历史销售
gen last_month_sales = sales_amount * (0.8 + runiform() * 0.4)
label variable last_month_sales "上月销售额"

gen sales_growth = (sales_amount - last_month_sales) / last_month_sales * 100
label variable sales_growth "销售增长率（%）"

*------------------------------------------------------------------------------
* 3. 探索性数据分析
*------------------------------------------------------------------------------

display ""
display "探索性数据分析"
display "=" * 80

* 描述性统计
summarize sales_amount unit_price product_quality market_share ///
    marketing_spend brand_strength

* 相关性分析
display ""
display "关键变量相关性："
correlate sales_amount unit_price product_quality market_share ///
    marketing_spend brand_strength

*------------------------------------------------------------------------------
* 4. 数据可视化
*------------------------------------------------------------------------------

display ""
display "生成探索性图表"
display "-" * 80

* 4.1 销售额分布
histogram sales_amount, normal ///
    title("销售额分布") ///
    xtitle("月销售额（千元）") ytitle("频数") ///
    scheme(s2color)
graph export "output/cases/figures/case03_sales_distribution.png", replace

* 4.2 价格与销售关系
scatter sales_amount unit_price, ///
    title("单价与销售额关系") ///
    xtitle("单价（千元）") ytitle("销售额（千元）") ///
    scheme(s2color)
graph export "output/cases/figures/case03_price_sales.png", replace

* 4.3 营销投入与销售关系
scatter sales_amount marketing_spend, ///
    title("营销投入与销售额关系") ///
    xtitle("营销费用（千元）") ytitle("销售额（千元）") ///
    scheme(s2color)
graph export "output/cases/figures/case03_marketing_sales.png", replace

* 4.4 市场份额与销售关系
scatter sales_amount market_share, ///
    title("市场份额与销售额关系") ///
    xtitle("市场份额（%）") ytitle("销售额（千元）") ///
    scheme(s2color)
graph export "output/cases/figures/case03_marketshare_sales.png", replace

* 4.5 不同价格档次的销售额
graph box sales_amount, over(price_level) ///
    title("不同价格档次的销售额分布") ///
    ytitle("销售额（千元）") ///
    scheme(s2color)
graph export "output/cases/figures/case03_price_level_sales.png", replace

* 4.6 季节性销售模式
graph bar (mean) sales_amount, over(season) ///
    title("季节性销售模式") ///
    ytitle("平均销售额（千元）") ///
    blabel(bar, format(%9.0f)) ///
    scheme(s2color)
graph export "output/cases/figures/case03_seasonal_sales.png", replace

*------------------------------------------------------------------------------
* 5. 数据分割
*------------------------------------------------------------------------------

display ""
display "数据分割（80% 训练集，20% 测试集）"
display "-" * 80

set seed 20251103
gen random = runiform()
gen train = (random <= 0.8)
gen test = (random > 0.8)

tab train
display "训练集: " _N if train
display "测试集: " _N if test

*------------------------------------------------------------------------------
* 6. 模型1：OLS 回归
*------------------------------------------------------------------------------

display ""
display "模型1：OLS 线性回归"
display "=" * 80

* 训练模型
regress sales_amount unit_price product_quality market_share ///
    marketing_spend brand_strength promotion_intensity last_month_sales ///
    if train

* 保存模型
estimates store ols_model

* 预测
predict sales_pred_ols

* 评估（测试集）
display ""
display "测试集评估："
display "-" * 80

gen sq_error_ols = (sales_amount - sales_pred_ols)^2 if test
egen mse_ols = mean(sq_error_ols)
gen rmse_ols = sqrt(mse_ols)

gen abs_error_ols = abs(sales_amount - sales_pred_ols) if test
egen mae_ols = mean(abs_error_ols)

correlate sales_amount sales_pred_ols if test
local r2_ols = r(rho)^2

summarize rmse_ols mae_ols
display "R²: " %6.4f `r2_ols'

*------------------------------------------------------------------------------
* 7. 模型2：Lasso 回归
*------------------------------------------------------------------------------

display ""
display "模型2：Lasso 回归"
display "=" * 80

* 训练 Lasso
lasso linear sales_amount unit_price product_quality product_size ///
    product_performance market_share competitor_count brand_strength ///
    marketing_spend promotion_intensity last_month_sales ///
    if train, selection(cv) rseed(20251103)

* 查看选择的特征
display ""
display "Lasso 选择的特征："
lassocoef

* 预测
predict sales_pred_lasso

* 评估
gen sq_error_lasso = (sales_amount - sales_pred_lasso)^2 if test
egen mse_lasso = mean(sq_error_lasso)
gen rmse_lasso = sqrt(mse_lasso)

gen abs_error_lasso = abs(sales_amount - sales_pred_lasso) if test
egen mae_lasso = mean(abs_error_lasso)

correlate sales_amount sales_pred_lasso if test
local r2_lasso = r(rho)^2

display ""
display "测试集评估："
summarize rmse_lasso mae_lasso
display "R²: " %6.4f `r2_lasso'

*------------------------------------------------------------------------------
* 8. 模型3：Elastic Net
*------------------------------------------------------------------------------

display ""
display "模型3：Elastic Net"
display "=" * 80

* 训练 Elastic Net
elasticnet linear sales_amount unit_price product_quality product_size ///
    product_performance market_share competitor_count brand_strength ///
    marketing_spend promotion_intensity last_month_sales ///
    if train, alpha(0.5) selection(cv) rseed(20251103)

* 预测
predict sales_pred_enet

* 评估
gen sq_error_enet = (sales_amount - sales_pred_enet)^2 if test
egen mse_enet = mean(sq_error_enet)
gen rmse_enet = sqrt(mse_enet)

gen abs_error_enet = abs(sales_amount - sales_pred_enet) if test
egen mae_enet = mean(abs_error_enet)

correlate sales_amount sales_pred_enet if test
local r2_enet = r(rho)^2

display ""
display "测试集评估："
summarize rmse_enet mae_enet
display "R²: " %6.4f `r2_enet'

*------------------------------------------------------------------------------
* 9. 模型对比
*------------------------------------------------------------------------------

display ""
display "模型性能对比"
display "=" * 80

* 创建对比表
quietly summarize rmse_ols if test
local rmse_ols_val = r(mean)
quietly summarize mae_ols if test
local mae_ols_val = r(mean)

quietly summarize rmse_lasso if test
local rmse_lasso_val = r(mean)
quietly summarize mae_lasso if test
local mae_lasso_val = r(mean)

quietly summarize rmse_enet if test
local rmse_enet_val = r(mean)
quietly summarize mae_enet if test
local mae_enet_val = r(mean)

matrix results = J(3, 3, .)
matrix rownames results = "OLS" "Lasso" "ElasticNet"
matrix colnames results = "RMSE" "MAE" "R²"

matrix results[1, 1] = `rmse_ols_val'
matrix results[1, 2] = `mae_ols_val'
matrix results[1, 3] = `r2_ols'

matrix results[2, 1] = `rmse_lasso_val'
matrix results[2, 2] = `mae_lasso_val'
matrix results[2, 3] = `r2_lasso'

matrix results[3, 1] = `rmse_enet_val'
matrix results[3, 2] = `mae_enet_val'
matrix results[3, 3] = `r2_enet'

display ""
matrix list results

*------------------------------------------------------------------------------
* 10. 预测结果可视化
*------------------------------------------------------------------------------

display ""
display "生成预测结果图表"
display "-" * 80

* 10.1 预测 vs 实际（OLS）
scatter sales_amount sales_pred_ols if test || lfit sales_amount sales_pred_ols if test, ///
    title("OLS模型：预测 vs 实际销售额") ///
    xtitle("实际销售额") ytitle("预测销售额") ///
    legend(label(1 "数据点") label(2 "拟合线")) ///
    scheme(s2color)
graph export "output/cases/figures/case03_ols_prediction.png", replace

* 10.2 三个模型对比
twoway (scatter sales_amount sales_pred_ols if test, mcolor(blue) msymbol(circle)) ///
       (scatter sales_amount sales_pred_lasso if test, mcolor(red) msymbol(triangle)) ///
       (scatter sales_amount sales_pred_enet if test, mcolor(green) msymbol(square)) ///
       (lfit sales_amount sales_amount if test, lcolor(black)), ///
       legend(label(1 "OLS") label(2 "Lasso") label(3 "Elastic Net") label(4 "完美预测")) ///
       title("模型预测对比") ///
       xtitle("实际销售额") ytitle("预测销售额") ///
       scheme(s2color)
graph export "output/cases/figures/case03_model_comparison.png", replace

* 10.3 残差分析
gen resid_ols = sales_amount - sales_pred_ols if test

histogram resid_ols, normal ///
    title("OLS模型残差分布") ///
    xtitle("残差") ytitle("频数") ///
    scheme(s2color)
graph export "output/cases/figures/case03_residuals.png", replace

* 10.4 预测误差对比
graph box sq_error_ols sq_error_lasso sq_error_enet if test, ///
    title("模型预测误差对比") ///
    ytitle("平方误差") ///
    legend(label(1 "OLS") label(2 "Lasso") label(3 "Elastic Net")) ///
    scheme(s2color)
graph export "output/cases/figures/case03_error_comparison.png", replace

*------------------------------------------------------------------------------
* 11. 业务洞察
*------------------------------------------------------------------------------

display ""
display "=" * 80
display "业务洞察与建议"
display "=" * 80
display ""

* 关键驱动因素（基于 OLS 模型）
estimates restore ols_model

display "销售额关键驱动因素："
display "-" * 80

* 使用边际效应
margins, dydx(*) post

* 保存结果
estimates store ols_margins

* 绘制系数图
coefplot ols_margins, drop(_cons) xline(0) ///
    title("销售驱动因素分析") ///
    xtitle("对销售额的影响") ///
    scheme(s2color)
graph export "output/cases/figures/case03_drivers.png", replace

*------------------------------------------------------------------------------
* 12. 保存结果
*------------------------------------------------------------------------------

* 保存预测结果
preserve
keep if test
keep sales_amount sales_pred_ols sales_pred_lasso sales_pred_enet ///
    unit_price market_share marketing_spend
export delimited "output/cases/case03_predictions.csv", replace
restore

* 保存完整数据
save "data/cases/sales_forecasting_results.dta", replace

display ""
display "=" * 80
display "销售预测分析完成！"
display "=" * 80
display ""
display "最佳模型: " cond(`r2_ols' > `r2_lasso', "OLS", "Lasso")
display ""
display "输出文件："
display "  - 预测结果: output/cases/case03_predictions.csv"
display "  - 图表: output/cases/figures/case03_*.png (10个图表)"
display "  - 日志: output/cases/case03_sales_forecasting.log"

log close

